/*
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

#include <src/math/Operators/D_x_Delta_S.h>
#include <src/algorithm/Evolution_Algorithm/Crank_Nicolson.h>
#include <src/math/CauchyProblems/CauchyProblem_traits.h>


#if !defined(_Cauchy_Burger_p_)
#define _Cauchy_Burger_p_

#define NUMBER_OF_PARTICLES 10

template <class data, std::size_t Dim = data::Dim>
struct EScheme_Burger_p_interface {
	OS_STATIC_CHECK(false); // only dim == 1 or dim == 2.	
};

template <class data>
struct EScheme_Burger_p_interface<data,1> : public EScheme<Burger_p_t<data>>
{
	typedef  EScheme										Base;
	typedef typename EScheme_Burger_p_interface				type;
	typedef Tadmor<data>									tadmor_type;
	typedef typename data::value_type						value_type;
	typedef typename tadmor_type::smart_ptr_function_type	smart_ptr_function_type;
	typedef typename data::base_value_type					base_value_type;

	template <class D>
	struct v : public integrable<D,D>
	{
		struct int_v : public Base
		{
			virtual result_type operator()(argument_type from) const {
				return (0.5)/(from);
			};
		};
		virtual result_type operator()(argument_type from) const {
			return  (-0.5) / pow(from,2) ;
		};
		virtual smart_ptr<Base> get_integral() const {
			return smart_ptr<Base>(new int_v);
		};
	};
	template <class D>
	struct g : public integrable<D,D>
	{
		struct int_g : public Base
		{
			virtual result_type operator()(argument_type from) const {
				return ( std::pow(-from,3./2.))/( std::sqrt(2.)/3.);
			};
		};
		virtual result_type operator()(argument_type from) const {
			return ( std::pow(-from,1./2.))/( std::sqrt(2.)/2.);
		};
		virtual smart_ptr<Base> get_integral() const {
			return smart_ptr<Base>(new int_g);
		};
	};
	type() : vstar_f_(smart_ptr_function_type(new v<value_type>)), gstar_f_(smart_ptr_function_type(new g<value_type>))
	{
		u_.resize(2,0);
	};
	virtual void vstar(indice_type n) 
	{
		vstar_[n] = vstar_f_( Op_(0.,u_[n]) , Op_(0.,u_[n+1]) );
	};
	virtual void gstar(indice_type n) 
	{
		gstar_[n] = gstar_f_(vstar_[n]>>1,vstar_[n])+Delta_p_((u_[n]+u_[n+1])/2.)*viscosity_;
	};
	virtual OS_double Estar(indice_type n) {
		v<data> U;
		return jmmath::sum( U(Op_(0.,u_[n])));
	};
	virtual OS_double CFL(indice_type n) {
		OS_double breakpoint1 = jmmath::norm_infty(gstar_[n]);
		return .5/ get_state_variable()->size(1)/breakpoint1;
	};
	virtual bool criteria(indice_type n) {
/*		OS_double breakpoint1 = Estar(n+1);
		OS_double breakpoint2 = Estar(n);
		OS_double breakpoint3 = (breakpoint1-breakpoint2)/(delta_t_*breakpoint2);
		return std::abs(breakpoint3) < 1e-4;*/
		
		OS_double breakpoint4 = jmmath::norm22(gstar_[n] - Op_S_k1);
		OS_double breakpoint5 = jmmath::norm22(gstar_[n]);
		return breakpoint4/breakpoint5 < 1e-4;
	};
	virtual void update(){
		Base::update();
		Delta_p_.set_state_variable(get_state_variable());
		viscosity_ = get_viscosity();
	}; 
	base_value_type get_viscosity() {
		base_value_type nb = 1./get_state_variable()->size();
		return nb;
	};
	void set_viscosity(base_value_type viscosity) {
		viscosity_ = viscosity;
	};
	base_value_type viscosity_;


	tadmor_type vstar_f_;
	tadmor_type gstar_f_;
	Delta_p<data>	Delta_p_;
};


/*! \brief Burger equation - Particle Method
    This class initialize the Cauchy problem for the one dimensional Burger equation
\f[ \label{burger_u}
    \partial_t u = \partial_x f(u); u(0,x) = u_0(x) > 0, \int_{\RR} u_0 = 1
\f]
    with \f$x \in \RR, t\ge 0\f$.
Since this equation owns many weak solutions, we consider an entropy / entropy flux function \f$(U(u),F(u))\f$ such that the solution verifies
weakly
\f[ \label{burger_e}
    \partial_t U(u) + \partial_x F(u) \le 0
\f]
 We provide a one dimensional particle method of this Burger equation.
 We recall that the particle method for the Burger equation consists in considering the equivalent equation
\f[ \label{burger_S}
    \partial_t S =  \tilde{f}\left(\partial_y S \right); S(0,y) = S_0(y)
\f]
with \f$y \in [0,1]\f$ and \f$\tilde{f}(s):=sf(\frac{1}{s})\f$, \f$S(t,y)\f$ being the inverse of the cumulative of \f$u(t,x)\f$, that we denote \f$S^{-1}(t,x)=\int_{-\infty}^xu(t,s)ds\f$.
The equation \f$(\ref{burger_S})\f$ is equivalent to the equations (\ref{burger_u})-(\ref{burger_e}) selecting a solution \f$S\f$ verifying the corresponding entropy / entropy flux inequality
\f[ \label{burger_eS}
    \partial_t \tilde{U}(\partial_y S) + \partial_y \tilde{F}\left( \partial_y S\right) \le 0
\f]
where we denoted
\f[
    \tilde{U}(s) = sU(\frac{1}{s});\tilde{F}(s) = F(\frac{1}{s})
\f]
Note that, calling \f$v(s)=\tilde{U}'(s)\f$ the entropy variable, and defining
\f$g\f$ through the relation
\f[
    v(s) = \tilde{U}'(s);g(v(s)) = \tilde{f}(s)
\f]
the particle equation \f$(\ref{burger_S})\f$ may be written equivalently as
\f[ \label{burger_S_e}
    \partial_t S =  g\left(v\left(\partial_y S\right) \right); S(0,y) = S_0(y)
\f]
We approximate the equation (\ref{burger_S_e}) is approximated
with the following dynamical system of \f$N\f$ particles \f$ \{S_i \}_{i=1,..,N} \f$ :
\f[ \label{sd_burger}
    \frac{d}{dt} S_i = g^*\left( v(\delta_{\bar{y}}S_i),v(\delta_{y}S_i)) \right); S_i(0) = S_i^0
\f]
 where \f$\{S_i \equiv S_i(t)\}_{i=1,..,N}\f$ are the particles location,
\f$\{S_i^0\}_{i=1,..,N}\f$ is the optimal quantizer of \f$u_0\f$ : i.e. \f$S_i^0 = S^0 (\frac{i+1/2}{N}), i=0,..,N-1\f$,
where \f$S^0\f$ is the inverse function of the cumulative of \f$u_0\f$, that we denote \f$(S^0)^{-1}(x) = \int_{-\infty}^x u_0(s)ds\f$.
Finally \f$g^*(v_0,v_1)\f$ is the Tadmor flux associated to \f$g\f$, that is
\f[
   g^*(v_0,v_1) = \frac{\int_{v_0}^{v_1}g(s)ds}{v_1-v_0}
\f]
Indeed, it can be proven easily that this semi discrete scheme verifies the semi discrete entropy condition
\f[
    \partial_t \tilde{U}(\partial_y S_i) + \delta_y \tilde{G}^*\left( v(\partial_{\bar{y}} S_i),v(\partial_y S_i)\right) = 0, i=0,..,N-1
\f]
where \f$G^*(v_0,v_1)\f$ is the (Tadmor) flux associated to \f$g^*(v_0,v_1)\f$.
The hamilton system (\ref{sd_burger}) is approximated directly by the non linear system of equations
\f[ \label{fd_burger}
    \delta_t S_i^n = g^*\left( v^{*,n+1/2}_{i-1,i},v^{*,n+1/2}_{i,i+1}) \right); S_i(0) = S_i^0
\f]
where
\f[
    v^{*,n+1/2}_{i,i+1} = v^*\left(\delta_{y}S_i^n,\delta_{y}S_i^{n+1} \right)
\f]
and \f$v^*(v_0,v_1)\f$ is the Tadmor flux associated to \f$v\f$, that is
\f[
   v^*(v_0,v_1) = \frac{\int_{v_0}^{v_1}v(s)ds}{v_1-v_0}
\f]
The choice of such a scheme is motivated by the fact that it verifies
the entropy evolution property
\f[
    \delta_t \sum_i \tilde{U}(\delta_{y}S_i^n) = 0
\f]
\begin{proof}
\f[
    \delta_t \tilde{U}(\delta_{y}S_i) = v^{*,n+1/2}_{i,i+1} \delta_t\delta_{y}S_i^n = v^{*,n+1/2}_{i,i+1} \delta_{y} g^*\left( v^{*,n+1/2}_{i-1,i},v^{*,n+1/2}_{i,i+1} \right)
\f]
Using the property of Tadmor fluxes, we have
\f[
    \delta_t \tilde{U}(\delta_{y}S_i) = \delta_{y}G (v^{*,n+1/2}_{i-1,i})
\f]
Summing over all indices \f$i\f$ will cancel the right hand term.
\end{proof}
Thus this scheme suits well to our purposes.
To solve the system \f$\ref{fd_burger}\f$, we use the Entropic Scheme fixed point algorithm : first initiate
\f[
    S_i^{0,n+1}=S_i^{n}.
\f]
then loop for \f$k \ge 1\f$:
\f[
    v^{*,k,n+1/2}_{i,i+1}:=v^*\left(\delta_{y}S_i^{n},\delta_{y} S_i^{k-1,n}\right);S_i^{k,n+1} = S_i^{n} + \tau^n g^*\left( v^{*,n+1/2}_{i-1,i},v^{*,n+1/2}_{i,i+1}) \right)
\f]
Until a convergence criteria holds, a generic and natural convergence criteria being
\f[
    \|\frac{\tilde{U}(\delta_{y}S_i^n)-\tilde{U}(\delta_{y}S_i^{k,n+1})}{\tau^n}\|_{\ell_1} \le \epsilon \| \tilde{U}(\delta_{y}S_i^n) \|_{\ell_1}
\f]

To illustrate numerically, we consider the classical Burger equation with a quadratic Entropy
\f[
    \partial_t u + \frac{1}{2}\partial_x u^2 = 0; \frac{1}{2} \partial_t u^2 + \frac{1}{3} \partial_x u^3 \le 0
\f]
This corresponds to the particle system
\f[
    \partial_t S + \tilde{f}(\partial_y S)=\partial_t S+\frac{1}{2\partial_y S}=0; \partial_t \frac{1}{2\partial_y S}  + \partial_y \frac{1}{3 (\partial_y S)^3} \le 0
\f]
This yields the expressions
\f[
    v(s) = -\frac{1}{2s^2}; g(s) = -\frac{\sqrt{2s}}{2}
\f]
and the Tadmor fluxes are
\f[
    v^*(v_0,v_1) = -\frac{1}{2v_0v_1}; g^*(v_0,v_1) = -\frac{\sqrt{2}}{3}\frac{v_0+v_1+\sqrt{v_0v_1}}{\sqrt{v_0}+\sqrt{v_1}}
\f]


Here, the initial condition correspond to a cauchy data given by a concentrated gaussian
\f$u_0(x) = N(\sigma,x) = \frac{1}{\sqrt{2 \sigma}} \exp {\frac{-x^2}{\sigma}}\f$, where \f$\sigma=0.1\f$.
This leads to initial condition given by
\f$\{S_i^0 = P^{-1}(\sigma,\frac{i+1/2}{N})\}_{i=1,..,N-1}\f$, where \f$P^{-1}(\sigma,x)\f$ is the inverse of the gaussian cumulative.
*/
		template <class data, std::size_t Dim = data::Dim, class Interface = Graphical_Interface>
		class Cauchy_Burger_p : public CauchyProblem_traits< EScheme_Burger_p_interface < data >, Interface>
		{
			OS_STATIC_CHECK(data::Dim == 1);

			public :
			virtual smart_ptr<interface_type> Create_Instance() {
				return smart_ptr<interface_type>(new Cauchy_Burger_p<data,Dim>);
			};

			Cauchy_Burger_p(){
				initial_conditions(NUMBER_OF_PARTICLES);
			};
			virtual void set_dim(OS_size dim_d)
			{
				OS_DYNAMIC_CHECK(dim_d == 1,"Cauchy_Burger_p : Only dimension 1");
				initial_conditions(NUMBER_OF_PARTICLES);
			};
			OS_size get_dim() {return 1;};
			void initial_conditions(OS_size nb_part)
			{
				set_initial_conditions(InverseCumulativeFunctor<data>(0.,1.)(nb_part));
			};
			virtual typename interface_type::p_array_2 get_representation() {
				return representation_2<result_type,result_type::Dim,Particle_Operator>()(get_solution());
			};
			virtual std::string get_name() {return "Burger Equation. Density Tracking.";};
			virtual void set_viscosity(OS_double viscosity)
			{
				get_algorithm().set_viscosity(viscosity);
			};
		};



#undef NUMBER_OF_PARTICLES
#undef DIM_D
#endif